***************************************************************************************************************
* Name: TableA6
* Description: runs biprobit model with exogenous variables and computes Vuong test and other stats reported in A6 (see output)
* Reference: "Competition and Civilian Victimization"
***************************************************************************************************************

use "violence_t_export.dta",clear
keep if victim_farc~=.&victim_paras~=.

local controls "nbi_t royalties_t coca_t share_left lpobl_tot_t i.time army_dist gini_i dmr evlp gcaribe gpacifica gorinoquia gamazonia"
biprobit (victim_farc = farc_dist `controls') (victim_paras= paras_dist `controls')

gen q_farc=1
replace q_farc=-1 if victim_farc==0
gen q_paras=1
replace q_paras=-1 if victim_paras==0
gen rhos=q_farc*q_paras*e(rho)
predict xb_farc,xb1
predict xb_paras,xb2
gen log_Lbi=ln(binormal(q_farc*xb_farc,q_paras*xb_paras,rhos))

drop xb_farc xb_paras q_paras q_farc

*****Log likelihood definition*****

global vF "victim_farc"
global vP "victim_paras" 

*Victimization un-restricted
capture program drop lfsec_ur
program lfsec_ur
  version 13.1
  args lnf xbF xbP
  tempvar pF pP
  quietly {
     gen double `pF' = 1/(1+exp(-`xbF'))
     gen double `pP' = 1/(1+exp(-`xbP'))
     replace `lnf' = $vP* ln(`pP')+(1-$vP)*ln(1-`pP') + $vF*ln(`pF')+(1-$vF)*ln(1-`pF') 
	}
end

*Two stage estimation
 capture program drop game2stage
 program game2stage, eclass
	syntax [if]
	marksample touse
	tempvar predF predP
	local controls "nbi_t royalties_t coca_t share_left lpobl_tot_t i.time army_dist gini_i dmr evlp gcaribe gpacifica gorinoquia gamazonia"
	
	quietly {
		npregress kernel victim_farc `controls' farc_dist if victim_paras~=., estimator(constant) noderivatives
		predict `predF'
		npregress kernel victim_paras `controls' paras_dist if victim_paras~=.,estimator(constant) noderivatives
		predict `predP'
	}
	
	gen pF=`predF'
	gen pP=`predP'
		
	corr  pF victim_farc pP victim_paras
	ml model lf lfsec_ur (eq1: victim_farc = nbi_t royalties_t coca_t share_left lpobl_tot_t army_dist  ///
						 dmr gini_i time evlp gcaribe gpacifica gorinoquia gamazonia farc_dist pP) (eq2: nbi_t royalties_t coca_t ///
						share_left lpobl_tot_t army_dist dmr gini_i time ///
						 evlp  gcaribe gpacifica gorinoquia gamazonia paras_dist pF)
	capture noisily ml maximize, diff 
	display "maximize" _rc

end

quietly: game2stage

predict xbF,equation(eq1)
predict xbP,equation(eq2)
gen double pF_aux = 1/(1+exp(-xbF))
gen double pP_aux = 1/(1+exp(-xbP))
gen log_LL = victim_paras* ln(pP_aux)+(1-victim_paras)*ln(1-pP_aux) + victim_farc*ln(pF_aux)+(1-victim_farc)*ln(1-pF_aux) 


**Voung and Clarck Test**

gen LR=log_LL-log_Lbi
gen LR_sq=LR^2
gen Ones=1


gen ind_d=0
replace ind_d=1 if LR>0

gen LR_corr=(log_LL-(34/(2*611)*ln(611)))-(log_Lbi-(33/(2*611)*ln(611)))

gen ind_dc=0
replace ind_dc=1 if LR_corr>0

mkmat LR,mat(lr)
mkmat LR_sq,mat(lr_sq)
mkmat Ones,mat(ones)
mkmat log_LL,mat(log_ll)
mkmat log_Lbi,mat(log_lbi)

matrix w_sq=inv(ones'*ones)*(lr_sq'*ones)-(inv(ones'*ones)*(lr'*ones))*(inv(ones'*ones)*(lr'*ones))
matrix lrr=lr'*ones
matrix N=ones'*ones
matrix log_LL_sum=log_ll'*ones
matrix log_Lbi_sum=log_lbi'*ones

scalar w_sq=w_sq[1,1]

*We correct by number of params
scalar lrr=lrr[1,1]
scalar N=N[1,1]
scalar log_LL_sum=log_LL_sum[1,1]
scalar log_Lbi_sum=log_Lbi_sum[1,1]

*p value Voung
di 1-normal((1/sqrt(N))*(lrr-((34/2)*ln(N)-(33/2)*ln(N)))/(w_sq^.5)) 

tab ind_dc
*p value clarke test
di 1-binomial(611,327,0.5)


*BIC
di -2*log_LL_sum+ln(N)*34
di -2*log_Lbi_sum+ln(N)*33

*AIC
di -2*log_LL_sum+2*34
di -2*log_Lbi_sum+2*33

*Checking Kurtosis
summarize LR ,detail
